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ABSTRACT 

Radiation dominated accretion disks are likely to be subject to the "photon bubble" instability, which 
may lead to strong density inhomogeneities on scales much shorter than the disk scale height. Such 
disks - and magnetized, radiation-dominated atmospheres in general - could radiate well above the 
Eddington limit without being disrupted. When density contrasts become large over distances of order 
the photon mean free path, radiative transfer cannot be described adequately using either the stan- 
dard diffusion approximation or existing prescriptions for flux-limited diffusion. Using analytical and 
Monte Carlo techniques, we consider the effects of strong density gradients deep within radiation- and 
scattering-dominated atmospheres. We find that radiation viscosity - i.e., the off-diagonal elements of 
the radiation stress tensor - has an important effect on radiative transfer under such conditions. We 
compare analytical and numerical results in the specific case of a plane-parallel density wave structure 
and calculate Eddington enhancement factors due to the porosity of the atmosphere. Our results can 
be applied to the study of dynamical coupling between radiation forces and density inhomogeneities in 
radiation dominated accretion disks in two or three dimensions. 

Subject headings: radiative transfer - accretion, accretion disks - stars: atmospheres - methods: 
analytical, numerical 



1. INTRODUCTION 

Radiation-dominated atmospheres of accretion disks 
and massive stars permeated by a moderately strong mag- 
netic field are susceptible to the "photon bubble" instabil- 
ity (Arons 1992; Gammie 1998; Begelman 2001). In this 
process, high density regions tend to be pulled downward 
along the field lines and low density regions are pushed 
upward by radiation forces. The gas elements accelerated 
by radiation forces enter density maxima, where radia- 
tion forces decrease, and then progress downward again 
completing the cycle. Because subsequent acceleration 
episodes are increasingly large and magnetic tension pre- 
vents high density regions from spreading sideways, the 
density contrast increases. This leads to density inho- 
mogeneities on scales shorter than the characteristic scale 
height of the accretion disk or stellar atmosphere. Un- 
der such conditions, radiation tends to bypass high den- 
sity regions and travel more freely through tenuous ones. 
If the low- and high-density regions are dynamically cou- 
pled and most of the mass is in the high density phase, 
then the flux necessary to support the atmosphere against 
gravity can exceed Eddington limit. Photon bubble in- 
stability may be applicable to objects suspected of hav- 
ing super-Eddington luminosities, such as "ultraluminous 
X-ray sources" (Begelman 2002) and narrow-line Seyfert 
1 galaxies (King and Puchnarewicz 2002). It has also 
been argued that the instability may occur in simple non- 
magnetized Thomson atmospheres as they approach the 
Eddington limit (Shaviv 2001b). The origin of this insta- 
bility and the statistical properties of the inhomogeneities, 
such as their size relative to the size of the atmosphere, 
are different than that predicted by the magnetic photon 



bubble instability. This mechanism has been invoked to 
explain the discrepancy between relatively high outburst 
luminosities and relatively low outflow velocities in novae 
(Shaviv 1998, 2000, 2001a). 

Radiative transfer in a medium where parts of the gas 
are optically thin, such as the upper regions of atmospheres 
in particular, cannot be treated using the standard diffu- 
sion approximation. The application of this approxima- 
tion in the optically thin regime may lead to the radi- 
ation propagation rate exceeding the free-streaming rate 
|F| = cu, where F and u are the radiation flux and energy 
density, respectively. Previous work on radiative trans- 
fer under such conditions focused on the development of 
flux- limited diffusion approximations (e.g., Levermore and 
Pomraning (1981); Mcha and Zylstra (1991); Anile and 
Romano (1992)). A serious limitation of these methods 
is that they are applicable only in cases where the angu- 
lar distribution of the specific intensity is a slowly varying 
function of space and time. 

In this paper, we focus on radiative transfer deep within 
atmospheres (i.e., at high optical depth), where the effects 
of flux- limited diffusion are less important. But we con- 
sider the case where large density gradients - on scales 
of order the photon mean free path - lead to rapid fluc- 
tuations in the angular distribution function. It is easy 
to see why the standard diffusion approximations (with 
or without flux-limiting) fail under these circumstances. 
The standard diffusion approach predicts that the flux re- 
sponds instantaneously to local changes in density. But 
in reality the flux can only respond to inhomogeneities 
provided that the density changes on scales much larger 
than the photon mean free path. When density inhomo- 
geneities are optically thin, radiation does not "see" any 
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density contrasts. Using a simple ansatz, we have derived 
a modified diffusion equation and show that the effects of 
"photon viscosity" (i.e., off-diagonal elements of the ra- 
diation stress tensor) play an important role in radiative 
transfer in this case. Using analytical and Monte Carlo 
techniques we find that our analytical approach is much 
more accurate than standard diffusion or multi-stream ap- 
proximations. Our results can be applied to the study of 
dynamical coupling between radiation forces and density 
inhomogcncitics in radiation-dominated accretion disks. 

The structure of this paper is as follows. In section 2 
we discuss three-dimensional radiative transfer and show 
that the "photon viscosity" terms are important. In sec- 
tion 3 we constrain ourselves to the case of periodic, pla- 
nar density waves embedded deep within an atmosphere 
and derive the corresponding Eddington enhancement fac- 
tors under the assumption of global dynamical equilib- 
rium. In section 4 we compare analytical formulas with 
Monte Carlo calculations and discuss our results. We pro- 
pose that our method could be incorporated into existing 
radiation hydrodynamics codes such as the RHD module 
for ZEUS developed by Turner and Stone (2001). 

2. RADIATIVE TRANSFER IN AN INHOMOGENEOUS 
ATMOSPHERE 

We consider radiative transfer within infinite highly in- 
homogeneous atmospheres where the effects of flux- limited 
diffusion arc very much reduced. Moreover, in a porous 
atmosphere, the radiation diffusion timescale is most 
likely much shorter than the characteristic gas dynami- 
cal timescale, i.e., radiation is probably not trapped by 
the motion of gaseous inhomogcncitics (Begelman 2001). 
Therefore, for the purpose of calculating radiative transfer, 
we can neglect the time-dependence and motion of the gas 
distribution. We approximate the intensity distribution as 



/(x,0) = /o(x,0) + — r!-F(x), 
47r 



(1) 



where F is the local flux vector, f2 is the directional unit 
vector, and x is the position vector. This approxima- 
tion means that the first nontrivial correction to inten- 
sity 7o(x, fi) is symmetric with respect to the direction 
defined by the local flux. This has to be contrasted with 
a familiar case of a plane-parallel atmosphere. In such 
a case, it is customary to assume that the directional 
distribution of intensity /(x, Cl) is symmetric relative to 
the vertical direction which coincides with the direction 
of the flux vector. However, in the case of highly inho- 
mogencoiis atmosphere, radiation bypasses denser regions 
and the flux vector can rapidly change its orientation and 
only the volume-averaged flux is vertical. Thus, in our 
approach, the local "symmetry axis" is changing direction 
throughout the atmosphere and no global symmetry is re- 
quired. We also assume that all the odd moments of /q 
vanish, i.e., / tlilod^l = J CliCljClklod^ = 0, but otherwise 
make no assumptions about the directional dependence of 
Iq- The equation of radiative transfer for a 3D scattering 
atmosphere reads 



-Cl-VI = -I + — 

a 47r 



Idfl, 



(2) 



where a = pK is the scattering coefficient. Taking zeroth 
and first moments of eq. (2) we obtain 



and 



V-F = 



F = V-T, 

a 



(3) 



(4) 



where T is the radiation stress tensor, the components of 
which are 

1 



T 



(5) 



The closure relation for in terms of Fi and J = 
J Idfl can be obtained by calculating the second mo- 
ment of eq. (2), assuming the form of the intensity given 
by eq. (1). This leads to the equation for the radiation 
stress tensor: 



«_ 1 fia + fi) for«=j 
3 5(7c \ dx3 ' ax^ 1 ^ 

for I 7^ j 



1 ( dF\ 



dx^ 



(6) 



where u = A'kJ / c is the energy density and J is the mean 
intensity. Note that (i) the diagonal terms Tu may be dif- 
ferent from one another, and (ii) the off-diagonal elements 
of the radiation stress tensor do not vanish. The first point 
implies that this approach incorporates variable Edding- 
ton factors. The off-diagonal elements are responsible for 
"photon viscosity" and are non-zero even though bulk gas 
motions in the atmosphere were assumed to be negligible. 
This is because the "photon fluid" is moving through the 
gas, and may exert shear stresses. Substituting eq. (6) 
into eq. (4) and using eq. (3), we obtain: 



Fi 



c du 
3(7 dx^ 
1 da 
a dx 
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dx^ dx* 



(7) 



dx^ 
2 da dFi' 
\ wju- Kjju J a dx* dxi 
In the above equation, summation is imposed only over 
repeated j indices. Equation (7) and equation (3) are the 
governing equations of radiative transfer in our approxi- 
mation. 

3. 2d density inhomogeneities 

In order to illustrate our method, we now consider the 
simplified case of a plane-parallel wave density pattern (see 
Fig. 1). The scattering coefficient is assumed to be only a 
function of the distance ^ perpendicular to the slabs, i.e., 
a = pK = cr(^) with ^ = iJ,x+{l—iJ^y/^z, where n = cosip. 
We also assume that d/dy = and that the components of 
flux depend only on ^. We consider an atmosphere which 
is in global (i.e., volume-averaged) hydrostatic equilibrium, 
where radiation pressure balances gravity —gz. This im- 
plies that the gradient of energy density must be of the 
form 



v« = «'(Ove 



3(a) 



(8) 



where a prime denotes differentiation with respect to ^, 
(it) is the volume average of the scattering coefficient, and 
■pEdd = 9c/k is the Eddington flux. From equation (7) and 
equation (3) we obtain: 
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Edd -I- — 2 — 1 — I ) (9) 

where Q;^(/i) = 10/(1 + 4/^^ - Integrating eq. (3) 

over ^ we obtain 

^iF, + (1 - ^2)1/2^^^ ^ ^^^^^ = (1 _ ^2)1/2^^^^ (^0) 

where Fq is the integration constant. Using eq. (10) to 
eliminate Fx from eq. (9) we get 



l_ fFl\' 

v2. 



= (1 - /x2)Fo + /x^^Jkdd + ^ ( — 1 • (11) 

Multiplying equation (11) by cr and then volume- averaging 
it and demanding that the solution be bounded, we have 

{<j)Fmd = W = (1 - m')(<t)Fo + M'(^)i^Edd, (12) 
where the first equality in equation (12) comes from the 
requirement that the atmosphere be in global hydrostatic 
equilibrium. Therefore, Fq = F^dd and the final equation 
for the vertical flux is 



F, 



(13) 



3.1. Eddington enhancement factor 

Using equation (13) we can now calculate the Eddington 
enhancement factor I = (Fz/FEdd)- We can simplify and 
non-dimensionalize equation (13) by defining 



^Edd 

Thus, the Eddington factor is 



(l^A 



(14) 



;^(i-m') + (/), (15) 

where (/) is the volume average of /. Defining optical 
depth dr = ad^, letting the prime now denote differen- 
tiation with respect to r, and normalizing the scattering 
coefficient to its mean value (i.e. a ^ a/ (a)), we obtain a 
simplified form of equation (13) 



-f" 2 f 

J -a .f 



(16) 



We now consider a periodic slab model in which the scat- 
tering coefficient is given by: 



(7 



• if — n < T < 0, region 1 
(72 if < T < T2, region 2 



(17) 



The solutions to equation (16) in regions 1 and 2 are 
/i,2(t) = ai,2 sinh(aT) + 61,2 cosh(ar) H , (18) 

0'1,2 

where ai.2 and 61.2 are the integration constants. Func- 
tions / and /' have to be continuous across each slab 
boundary. Therefore the matching conditions are: 

/i(0) = /2(0), /{(0) = /^(0) (19) 
/i(-ri) = /2(t2), /{(-Ti) = /^(r2). (20) 



Using the above matching conditions to derive the inte- 
gration constants 01^2 and &i^2 and then volume averaging 
the solution for /, we obtain the expression for (/) 



6+6 



(Tl (72 



2 (Acr)^ sinha;sinhy 
sinh(a; -|- y) 



'1"2 



(21) 



where x = aTi/2, y = ar^ft, Act = ct2 — CTi, and 
6,2 = Ti,2/o'i,2- The Eddington enhancement factor I fol- 
lows directly from equations (21) and (15). 

4. COMPARISON WITH MONTE CARLO SIMULATIONS 

We now compare our analytical results with Monte 
Carlo simulations. The setup of the numerical experi- 
ment was as follows. We instantaneously injected a large 
number of photons in the equatorial plane (i.e., 2; = 0) 
of a very flat, three-dimensional box (i.e., z € [— 2:0,^0]; 

G [—Wo, Wo], where Zq <^ Wq). Densities and height 
of the computational box were chosen in such a way as 
to assure that the optical thickness in the vertical direc- 
tion would always be large throughout the box. The ini- 
tial photon angular distribution was uniform. Although 
we included the anisotropy due to the Thomson scatter- 
ing cross section, we found this effect to have no influ- 
ence on our final results. We followed the trajectories of 
all photons and calculated photon travel times between 
scatterings. Momentum transfer for every scattering was 
calculated using the method of weights (Pozdnyakov at 
al. 1983). We then computed the force exerted on the 
atmosphere as a function of the time delay following the 
instantaneous photon injection. Of course, as photons dif- 
fuse out of the atmosphere, the force exerted on the gas 
gradually declines. Therefore, the total force, correspond- 
ing to a continuous photon flux, was calculated by super- 
posing many such time-dependent force distributions due 
to groups of photons injected (instantaneously) at uniform 
time intervals. The total force exerted on the atmosphere 
was characterized by a gradual increase with time followed 
by a flat maximum. The Eddington enhancement factor is 
then given by the ratio of the "saturated" total force (i.e., 
total, constant force at late times) exerted on a homoge- 
neous atmosphere, to the total force acting on an inhomo- 
geneous atmosphere characterized by the same mean den- 
sity. Note that the Eddington enhancement factor, defined 
in this way, can also be interpreted as the ratio of fluxes 
necessary to exert the same amount of force on an inhomo- 
geneous atmosphere as on the corresponding homogeneous 
atmosphere. This ratio is the same even if the homoge- 
neous atmosphere is sub-Eddington, i.e., is only partially 
supported against gravity by radiation. 

Figure 2 shows the Eddington enhancement factor for 
variable inclination of slabs (upper panels) and for chang- 
ing density contrast of vertical slabs (lower panels). In all 
cases, the Thomson depth of the high density slabs is 
constant but the optical depth across low density regions 
increases from values < 1 to Tji > 1 from left to right 
(sec caption of Fig. 2 for more details). As expected, the 
Eddington factor increases as the slabs rotate toward the 
vertical direction because the atmosphere effectively be- 
comes more porous (see upper panels). When the slabs 
are vertical, the flux enhancement factor increases as the 
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density contrast Uh/cfi becomes larger for constant mean 
density. This is due to the fact that the volume filling 
factor of the high-density gas decreases while that of the 
low-density gas increases, but the respective masses of the 
two density phases remain the same. Therefore, the mean, 
volume- weighted, flux is 

{F) = {l-f,)Fi+UFh^Fu (22) 

where fy is the volume filling factor of the dense gas and 
Fi and Fh arc the fluxes propagating through tenuous and 
dense regions, respectively. As the density contrast in- 
creases and decreases, radiation tends to "flow" pri- 
marily through the low density channels and, therefore, 
more flux is necessary to exert the same total force as in 
the homogeneous case because radiation interacts less ef- 
flciently with tenuous gas. Quantitatively, in the diffusion 
limit, we have (Shaviv 1998) 

\p/ \<yi Oh) ii +6i 

When most volume is in the low-density phase but most 
mass is in the high-density phase (i.e., ^ or ^/j/^; ^ 
0), then I ?» fvPh/pi if fv > Pi/Ph- This qualitatively 
explains why I decreases with > 1 at constant density 
contrast (Jh/oi (cf. third and fourth columns on Fig. 2). 
At small optical depth ti , equation (23) would lead to very 
inaccurate answers. For example, equation (23) predicts 
I ~ 23 for vertical slabs with t; = 0.1 and at/cri = 100, 
compared to the actual value I ^ 5 and our analytic result 
^ 4 from eq. (21) (cf. lower left panel). This discrep- 
ancy is due largely to neglect of the anisotropy of the 
radiation field, whereas our approach gives much more 
accurate results even in such an extreme case. More- 
over, note that the "anisotropy term" in our expression 
for the flux enhancement factor, which is proportional to 
(Acr)^ = {(Th — oi)^, vanishes for large Thomson depths 
and thus equations (21) and (15) reduce to equation (23) 
in the diffusion limit. We also considered "multi-stream" 
approximation schemes in order to account for the radi- 
ation anisotropy, but found the "intensity moment" ap- 
proach developed here to be in significantly better agree- 
ment with Monte Carlo simulations. 



as the RHD module for ZEUS (Turner and Stone 2001). 
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5. SUMMARY 

We have considered radiative transfer deep within ex- 
tremely inhomogeneous atmospheres, and have demon- 
strated that, under such conditions, radiation viscosity 
- i.e., the off-diagonal elements of the radiation stress 
tensor - plays an important role. Our approach is sig- 
nificantly more accurate than approaches based on the 
diffusion equation and multi-stream approximation. The 
technique developed here can be applied to the nonlin- 
ear evolution of radiation-driven instabilities in accretion 
disks. In particular, it can be used to study the dynamical 
coupling of matter and radiation in order to determine 
the characteristic length scales and density contrasts aris- 
ing from "photon bubble" instability. This, in turn, will 
permit a self-consistent determination of the magnitude of 
the Eddington enhancement factor in radiation-dominated 
accretion disks. We also suggest that our method could 
be incorporated into radiation hydrodynamics codes such 
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Fig. 1. — Density structure of an inhomogeneous atmosphere. Shaded and unshaded zones denote higher and lower density regions, 
respectively. The region shown is much smaller than the overall radiation pressure scale height I I . 
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Fig. 2. — Eddington enhancement factors L/i^Edd for = 20 and r; = 0.10, 0.33, 1.00, 3.00 (from left to right column, respectively). 
Upper panels show L/Lgdd ^ ^ function of /i = cos^ for scattering coefficients of low and high density regions equal to ct; = 10 and 
(Tfi = 1000, respectively. Lower panels show L/L^^^ as a function of Ohjai for vertical slabs. Dashed lines show analytical predictions from 
equations (15) and (21), while solid lines show results from Monte Carlo simulations. 



